Anisotropic Migdal-Eliashberg theory using Wannier functions 
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We combine the fully anisotropic Migdal-Eliashberg theory with electron-phonon interpolation 
based on maximally-localized Wannier functions, in order to perform reliable and highly accurate 
calculations of the anisotropic temperature-dependent superconducting gap and critical temperature 
of conventional superconductors. Compared with the widely used McMillan approximation, our 
methodology yields a more comprehensive and detailed description of superconducting properties, 
and is especially relevant for the study of layered or low-dimensional systems as well as systems 
with complex Fermi surfaces. In order to validate our method we perform calculations on two 
prototypical superconductors, Pb and MgB2, and obtain good agreement with previous studies. 



I. INTRODUCTION 

The prediction of superconducting properties such as 
the critical temperature or the superconducting energy 
gap remains one of the outstanding challenges in mod- 
ern electronic structure theory.— ~— Owing to the com- 
plex nature of the superconducting state, a quantita- 
tive understanding of the pairing mechanism in conven- 
tional superconductors requires a very detailed knowl- 
edge of the electronic structure, the phonon dispersions, 
and the interaction between electrons and phonons. In 
conventional superconductors below the critical temper- 
ature electron pairing results from a subtle interplay be- 
tween the repulsive Coulomb interaction and the attrac- 
tive electron-phonon interaction. Starting from the sem- 
inal work of Bardeen, Cooper, and Schricffer (BCS) 8 
several approaches to the calculation of the supercon- 
ducting properties have been proposed, ranging from 
semi-empirical methods such as the McMillan formula^ 
to first-principles Green's function methods such as the 
Migdal-Eliashberg (ME) formalism ) 10 ! 11 and more re- 
cently also methods based on the density-functional the- 
ory concept, such as the density functional theory for 
superconductors (SCDFT)^ 

The vast majority of current investigations rely on the 
semi-empirical McMillan's approach^ In this approach 
the entire physics of the electron-phonon interaction is 
condensed into a single parameter called the electron- 
phonon coupling strength A. The McMillan's method 
works reasonably well for conventional bulk metals and 
for anisotropic superconductors where the Fermi surface 
anisotropy is smeared out by impurities.—"^ However, for 
layered systems, systems of reduced dimensionality, and 
those with complex multi-sheet Fermi surfaces, a care- 
ful description of the pairing interactions is crucial and a 
proper treatment of the anisotropic electron-phonon in- 
teraction is required. The necessity of an anisotropic the- 
ory has clearly been demonstrated in two important cases 
such as magnesium diboridc, MgB2, and the graphite in- 
tercalation compound CaCg, using either the ME formal- 
ism or the SCDFT.i£-i£ 

Unfortunately, the lack of adequate computational 
tools prevents the research community from systemat- 
ically exploring the importance of anisotropy in exist- 



ing and well characterized superconductors, and for val- 
idating computational predictions based on the semi- 
empirical McMillan equation*^— This latter aspect is 
particularly relevant in view of the increasingly im- 
portant role that high-throughput materials design ap- 
proaches are acquiring in the community 23 ' 24 

The momentum-resolved superconducting gap and the 
quasiparticlc density of states near the Fermi surface 
can now be measured with unprecedented accuracy us- 
ing high-resolution angle-resolved photoemission spec- 
troscopy^ and scanning tunneling spectroscopy 2 ^ experi- 
ments. In this context, the anisotropic Migdal-Eliashberg 
formalism promises to be particularly useful for perform- 
ing direct comparisons between theory and experiment, 
and for helping establish unambiguously the symmetry 
of the order parameter. 

One critical point which arises when attempting to 
solve the Eliashberg equations of the ME theory or the 
Bogoliubov-de Gennes equations of the SCDFT is that 
both sets of equations suffer from a strong sensitivity to 
the sampling of the electron-phonon scattering processes 
in the vicinity of the Fermi surfaced The practical con- 
sequence of this sensitivity is that, in order to achieve 
numerical convergence, the electron-phonon matrix ele- 
ments must be evaluated for extremely dense electron 
and phonon meshes in the Brillouin zone. This long- 
standing difficulty was overcome in Ref. by developing 
an efficient first-principles interpolation technique based 
on maximally- localized Wannier functions (MLWF)<2£ In 
this method one takes advantage of the spatial local- 
ization of both electron and phonon Wannier functions 
in order to evaluate only a small number of electron- 
phonon matrix elements in the Wannier representation. 
These matrix elements are subsequently interpolated to 
arbitrary electron and phonon wavevectors in the Bloch 
representation using a generalized Fourier transform^ 
This method carries general validity and has been demon- 
strated in several other areas, from Fermi surface calcu- 
lations^ to the anomalous Hall effect^ and more re- 
cently for GW calculations^ A detailed introduction on 
Wannier-based interpolation methods can be found in 
Ref. [28|. The scheme of Ref. [27] is adopted in the present 
work since it provides a robust and efficient framework for 
developing an algorithm to solve the anisotropic Eliash- 
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berg equations. Our current implementation enables the 
calculation of the momentum- and band-resolved super- 
conducting gap using a very fine Brillouin zone sampling. 
Without our Wannier-based electron-phonon interpola- 
tion this operation would not be possible owing to the 
prohibitive computational cost. 

This manuscript is organized as follows. In Sec. [TTI wc 
review the Migdal-Eliashbcrg theory of superconductiv- 
ity. A description of the computational techniques under- 
pinning the electron-phonon interpolation implemented 
in the EPW code 3 -^ is given in Sec. IIIII In Sec. IIVI we re- 
port the numerical solutions of the Eliashberg equations 
for two prototypical superconductors, Pb and MgP>2. Fi- 
nally we present our conclusions and outlook in Sec. [Vj 



II. MIGDAL-ELIASHBERG FORMALISM 

A. General theory 

A quantitative theory of the superconducting energy 
gap can conveniently be formulated within the framework 
of the Nambu-Gor'kov formalism! 33 - 34 In this formalism 
one introduces a two-component field operator: 



(1) 



The component c k -|- (c_ k i ) of the operator destroys (cre- 
ates) an electron in a Bloch state of combined band and 
momentum index k (— k) and spin up (down). A general- 
ized 2x2 matrix Green's functions G is then introduced 
in order to describe electron quasiparticles and Cooper's 
pairs on an equal footing)^ 



G(k,r) = -<T T * k (r 



'* k (o)>, 



(2) 



where T T is Wick's time-ordering operator for the imag- 
inary time r and ^^(t) is obtained from Eq. (TT]) using 
the Heisenberg picture. The braces indicate a grand- 
canonical thermodynamic average. By replacing Eq. (TT]) 
inside Eq. © we find: 



G(k,r) = - 



(T T c kt (r) Ckt (0)) (T TCkt (r)c_ ki (0)) 
(T rC t (t)c_^(0)> 



(7VcL 14 (T)c kt (0)) 



(3) 

Here the diagonal elements correspond to the standard 
Green's functions for electron quasiparticles and describe 
the dynamics of single-particle electronic excitations in 
the material. On the other hand, the off-diagonal ele- 
ments represent Gor'kov's anomalous Green's functions 
F(k, r) and F*(k, r). These functions describe the dy- 
namics of Cooper's pairs and are related to the supercon- 
ducting energy gap.— The off-diagonal elements of the 
generalized Green's function in Eq. (j3J) become nonzero 
only below the critical temperature T c , marking the tran- 
sition to the superconducting state. 



The generalized Green's function G(k, r) is periodic 
in imaginary time, therefore it can be expanded using a 
Fourier series as follows: 



G(k,r) =T^e- i; "" T G(k, lW „), 



(4) 



where itu n = i(2n + 1)ttT (n integer) stands for the 
fermion Matsubara frequencies, and T is the absolute 
temperature. We use atomic units throughout the 
manuscript, therefore we set h = Hb = 1- Following 
Eq. (TJ| the matrix elements of the generalized Green's 
function read: 



G(k,iw„) 



G(k, iuj n ) F(k,iw n ) 
F*(k,ioj n ) -G(-k,-iw n ) 



(5) 



The study of the superconducting state involves the de- 
termination of the matrix Green's function in Eq. ([5]). 
This can be achieved using Dyson's equation: 



G 1 (k,iu> n ) = G iuj n ) - £(k,iu> n ), 



(0) 



where Go(k, iuj n ) is the electron Green's function for 
the normal state and S(k, iuj n ) is the self-energy asso- 
ciated with the pairing interaction. The normal-state 
Green's function is calculated by using the Kohn-Sham 
states from density-functional theory to represent single- 
particle excitations. If we denote by e k the Kohn-Sham 
eigenvalues measured with respect to the chemical po- 
tential, and we introduce the Pauli matrices: 



to = 



r-2 



1 

1 

- 
i 



n = 



T3 = 



1 

1 o 

l o 

o -1 



(7) 



then the normal-state matrix Green's function acquires 
the familiar form: 



G 1 (k, iu n ) = iuj„T - e k f 3 . 



(8) 



Within the Migdal-Eliashberg approximation the elec- 
tron self-energy leading to the superconducting pairing 
consists of two terms, an electron-phonon contribution 
S cp (k, iu) n ) an d a Coulomb contribution £ c (k, iLu n )& 



£(k,iu> n ) = S cp (k, iu> n ) + £ c (k,wj„), 



(9) 



with 



E ep (k, iuj n ) = - T f 3 G (k' , iu n >) f 3 

k'n' 

| D v (k — k', iuj 



lUJ r , 



(10) 



and 



£ c (Kioj n ) = -T]Tf 3 G od (kV Wn 0Wk-k'). (11) 
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In Eq. (HDJ A,(q,iw„) = 2w q „/ [(iw„) 2 - w 2 „] ig thc 
dressed propagator for phonons with momentum q and 
branch index u, and is the screened electron-phonon 
matrix element for the scattering between the electronic 
states k and k' through a phonon with wavevector q = 
k' — k, frequency w q „ and branch index v. In Eq. (jlip 
the V(k — k')'s represent thc matrix elements of the 
static screened Coulomb interaction between the elec- 
tronic states k and k'. 

In writing the electron self-energy of Eqs. (p?))- (fTT]) thc 
following approximations are made: (i) only the first 
term in the Feynman's expansion of the self-energy in 
terms of electron-phonon diagrams is included. This ap- 
proximation corresponds to Migdal's theorem^ and is 
based on the observation that the neglected terms are 
of the order (mc/M) 1 / 2 (m e being the electron mass 
and M a characteristic nuclear mass), (ii) Only the off- 
diagonal contributions of the Coulomb self-energy are re- 
tained. This is done in order to avoid double counting the 
Coulomb effects that are already included in Go(k, ioj n )£ 
(iii) The self-energy is assumed to be diagonal in the 
electron band index^ This should constitute a reason- 
able approximation for non-degenerate bands since the 
energy involved in thc superconducting pairing is very 
small, therefore band mixing is not expected. 

In the literature on the theory of superconductivity it 
is common practice to decompose the matrix self-energy 
S(k, iu>„) in a linear combination of Pauli matrices with 
three scalar functions as coefficients. Thc scalar functions 
are the mass renormalization function Z(k, iuj n ), the en- 
ergy shift x(k, iu) n ), and the order parameter </>(k, iu) n ), 
and the decomposition reads: 

S(k, iui n ) = iu>„ [1 - Z(k, iu n )] f + x(k, iw n )f 3 

+ ^(k,iw n )fi. (12) 

We note that here we have chosen the gauge where the 
order parameter <f> is set to zero.— By replacing Eqs. (JSJ, 
(fT2f inside Eq. (|6|) and solving for the matrix Green's 
function we obtain: 

G(k, iui n ) = - {iuj n Z(k, iw„)f + [e k + x( k , iu> n )] h 

+ 0(k,icj„)fi}/6(k,iw„), (13) 

where the denominator is defined as: 

0(k,iw n ) = [w n Z(k, iuj n )} 2 + [e k + x(k, iw„)] 2 

+ [<Kk,iu n )]*. (14) 

At this point the strategy is to impose self-consistency by 
replacing the explicit expression for the Green's function 
in Eq. (|T3|) inside the self-energy expressions Eqs. ([9])- 
(JXTJ) . After equating the scalar coefficients of the Pauli 
matrices this replacement leads finally to the anisotropic 
Eliashberg equations: 

Z(k,l(J n ) = 1 + —^- 2^ a(v a, . \ A ( k ' k > n ~ n h 



X(k, iuj n 



u n N F f-J 9(k',iuv) 

k'rr 



T ^ e k / + x(k',iuj n >) 
N F ^ ®{k',iu) n ,) 



A(k, k', n - n ) 



(16) 



Hk^ n ) = ^E@^r— ; y 

K'rv 

x[\{k,k',n-n') - N F V(k~k')}. (17) 

In Eqs. (fT"5")) - p7[) N F represents the density of states per 
spin at the Fermi level, and A(k, k',n — n') is an auxil- 
iary function describing thc anisotropic electron-phonon 
coupling and defined as follows: 



A(k,k', n-n') = [ 
Jo 



2co 

doj- 



(uj n - LU n >) 2 + LO : 



r a 2 F(k,k», 

(18) 

with a 2 F(k, k',uj) the Eliashberg electron-phonon spec- 
tral function: 

a 2 F(k, k', u) = A F l-9kkv| 2 <5(^ - Wk-*,„). (19) 



The superconducting gap is defined in terms of the renor- 
malization function and thc order parameter as: 



A(k, ioj n ) 



(k, lLO ri 



Z(k,iw n ) 



(20) 



(15) 



From Eqs. (fTT|) . (|20[) we see that the Eliashberg equations 
admit the trivial solution A(k, iu) n ) = at all tempera- 
tures. Thc highest temperature for which the Eliashberg 
equations admit nontrivial solutions A(k,iu n ) ^ de- 
fines the critical temperature T c . 



B. Standard approximations 

After having presented a concise derivation of Eliash- 
berg's equations in the previous Section, we now discuss 
technical aspects which need to be addressed in order to 
actually solve these equations. 

Since thc superconducting pairing occurs mainly 
within an energy window of width w p h around the Fermi 
surface (w p h being a characteristic phonon energy), it 
is standard practice to simplify thc Eliashberg equa- 
tions by restricting the description to electron bands 
near the Fermi energy.— ~—i2£ This simplification can be 
achieved in the formalism by introducing thc identity 
de'<5(e k / — e') = 1 on thc right hand side in Eqs. (fT5")) - 
(fT"7l) . Thc rapid changes of 0(k', u> n >) and the numerator 
of Eq. (fH)| with the energy e' can be integrated analyt- 
ically, while for the other quantities we can set e' to the 
Fermi energy since the associated variations take place 
on a much larger energy scale^£ Under this approxima- 
tion thc energy shift becomes x(k, ioj n ) = and only two 
equations are left to solve, one for the renormalization 
function and one for the order parameter (or equivalcntly 
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the superconducting gap): 

Z(k, iu> n ) = l + — y2 Wv " n '. A(k, k', n-n'), 
1,1- ^ R{k!,iu n ,) 



(21) 



Z(k,iw„)A(k,iw„) = irT^W k 



A(k',iay) 
R(k',iu} n >) 



x[A(k > k',n-n / )-JVpV(k-k / )], (22) 



where i?(k, iw n ) and Wk are given by: 



J2(k, iw„) = + A 2 (k,iw n ) and W k = <5(e k )/iV F . 

(23) 

Equations (f2~Tj) and (f2"2"|) form a coupled nonlinear system 
and need to be solved self-consistently at each tempera- 
ture T. The approximations leading to Eqs. (f2~Tj) and (|22|) 
imply that Z(k, iw n ) and A(k, iw n ) are only meaningful 
for the momentum/band index k at or near the Fermi 
surface. Away from the Fermi surface the energy depen- 
dence of these quantities is weak and is neg lected^ In 
addition Eqs. (f2"T j) - (|2"2"j) implicitly assume that the elec- 
tron density of states is approximately constant near the 
Fermi energy. This simplification may break down for 
materials with narrow bands or critical points in proxim- 
ity of the Fermi level 

In order to solve Eqs. (|2"TT) . (|22|) numerically it is nec- 
essary to truncate the sum over Matsubara frequencies. 
It is standard practice to restrict the sum to frequencies 
smaller than a given cutoff lo c , with the cutoff of the or- 
der of 1 eV and typically set to 4-10 times the largest 
phonon energy. In addition, it is convenient to introduce 
a dimcnsionless Coulomb interaction parameter /i* de- 
fined as the double Fermi surface (FS) average over k 
and k' of the term V(k - k') in Eq. ([22]) : 



,i c = iV F «V(k-k')» 



FS ■ 



(24) 



By performing the energy integral analytically up to the 
cutoff frequency it can be shown that the N-pV(k — k') 
term in Eq. (|2"2"j) can be replaced by the Morel- Anderson 
pseudopotential /x* given by:— 



/'c 



/'c 



1 + /Lt c \n(e F /u c ) 



(25) 



Following this replacement, /i* is used as a semi-empirical 
parameter in the subsequent numerical solution of the 
Eliashberg equations. For a large class of superconduc- 
tors Eqs. ([24" |) - (f25]l yields values of /j* in the range 0.1- 
0.2. However, it is clear by now that in several cases 
values of the Coulomb parameter outside of this range 
are necessary for explaining experimental data . 27 ' 39 ' 40 In 
addition, the anisotropic nature of the Coulomb interac- 
tion cannot be neglected for an accurate description of 
the superconducting properties . 12 ' 14 ' 15 ' 38 These observa- 
tions should make it clear that the simplification provided 
by Eqs. (12"4"1) - ([2"5"1) is not optimal, and a fully ab-initio 
approach to the solution of the Eliashberg equations is 



highly desirable A description of the electron-phonon 
and the electron-electron interactions on the same footing 
is achieved in the SCDFT approach^ 7 - While it is clear 
that the Eliashberg approach considered in this work 
should be extended in order to incorporate Coulomb ef- 
fects from first-principles, this is beyond the scope of the 
present investigation. 



C. Superconducting gap along the real energy axis 

In Eqs. (|13[) - (|2"3"|) the dynamical aspects of the super- 
conducting pairing are described using the imaginary 
Matsubara frequencies ius n . The reason for this choice 
is that the resulting formulation is computationally ef- 
ficient since it only involves sums over a finite number 
of frequencies. While the imaginary axis formulation is 
adequate for calculating the critical temperature as de- 
scribed in Sec. Ill A[ the calculation of spectral properties 
such as the quasiparticlc density of states requires the 
knowledge of the superconducting gap along the real fre- 
quency axis. 

It is in principle possible to calculate the supercon- 
ducting gap along the real axis, however this procedure 
involves the evaluation of many principal value integrals 
and hence is numerically demanding In this work we 
prefer instead to determine the solutions of the Eliash- 
berg equations on the real axis by analytic continuation 
of our calculated solutions along the imaginary axis. The 
analytic continuation can be performed either by using 
Pade approximants as in Refs. 42.43 or by means of an 
iterative procedure as in Ref. 0. 

The continuation based on Pade approximants involves 
a very light computational workload, however it is very 
sensitive to the numerical precision of the solutions on the 
imaginary axisi 42 i 43 ' 45 As a rule of thumb the analytic 
continuation based on Pade approximants exhibits the 
correct gross structure of the superconducting gap on the 
real frequency axis, however fine spectral features are not 
always captured completely. 

The iterative analytic continuation, on the other hand, 
is generally rather accurate but involves a high compu- 
tational workload. In fact, as shown in Ref. |4J, the iter- 
ative analytic continuation requires solving the following 
equations sclf-consistcntly: 



Z(k,w) = 1+t— yiWu " n [ , A(k,k>-mv) 

k'n' 

+ i- / du'T(u,u') ")TWva 2 F (k,k',u/) 



(u - uj')Z(k',cu - w') 
P(k',w-w') ' 



(26) 
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A(kVay) . 
X P(k',w) 



Z(k, w)A(k, w) = ttT ^ W k , [A(k, k', w - iw„>) - 

k'n' 

/CO 
dwT(w, w') ^ W k 'a 2 F(k, k', w') 
-oo k , 

P(k',w-a/) ' 1 j 

where the following quantities have been introduced: 

, w) [w 2 + A 2 (k, w)], (28) 

.. If LJ — U. 

r(w,w ) = - I tanh ~ 2T~ 



A(k, k',w - iw„) 



coth-1, (29) 
^F(k,k>Q ; (3Q) 



Ul — lbJ n 



2 F(k,k',-w) = -a 2 F(k,k',w). 



(31) 



In the case where the square-root on the right hand side 
of Eq. (|28jl is complex, the root with positive imaginary 
part is chosen. 

Once determined the mass renormalization function 
Z(k, lu) and the superconducting gap A(k, uS) on the real 
frequency axis, one can examine the poles of the diagonal 
component of the single-particle Green's function^ 



Gu(k, w ) 



u>Z(k,uj) + e k 



[ojZ^w)} 2 - e 2 -[Z(k,w)A(k,w)]2 



in order to obtain the quasiparticle energy E k : 

2 



El = 



Z(k,E k ) 



A 2 (k,£ k ). 



,(32) 



(33) 



At the Fermi level e k = and the quasiparticle shift is 
Ek = RcA(k, Ek). As a result the leading edge A k of the 
superconducting gap for the combined band/momentum 
index k at the Fermi surface is given by: 



Re[A(k,A k )] = A k . 



D. Isotropic approximation 



(34) 



For conventional bulk metals or superconductors where 
the Fermi surface anisotropy is either weak or smeared 
out by impurities, it is possible to resort to a simplified 
isotropic formulation of the Eliashberg equations. Such 
formulation is obtained from Eqs. (f2Tj) . (f22j) by averaging 
k over the Fermi surface. We obtain: 



Z(iuj r , 



1 



ttT ■ 



0J n , R{lUJn') 



A(n-n'), (35) 



Z(iu n )A(iuj n ) = ttT V [A(n-n') - Mc *] , (36) 

TV 

where R(iuj n ) and X(n — n') are given by: 



R(iu n ) = v^ 2 + A 2 (iw„), (37) 



\(n-n')= du- ; „ , 38 

Jo K ~ w «') 2 + ^ 

and a 2 F(uj) is the isotropic Eliashberg spectral function: 

a 2 F(u) = Wk^a 2 F(k,k». (39) 
k , k ' 

The isotropic Eliashberg equations on the real axis can be 
obtained similarly by starting from Eqs. (f2l)]) . (f2"T)) . From 
the isotropic superconducting gap on the real axis we can 
obtain the normalized quasiparticle density of states in 
the superconducting state Ng(u)): 



N P 



Re 



y/0J 2 - A 2 (u) 



(40) 



III. COMPUTATIONAL METHODOLOGY 
A. Electron-phonon Wannier interpolation 

We now describe the combination of the anisotropic 
Eliashberg formalism of Sec. Ill Al with the electron- 
phonon Wannier interpolation of Ref. [I?]. The numer- 
ical solution of Eqs. ([2"T j) .([2"2" ]) and Eqs. ([2"g j) .([2"7 ]) requires 
an extremely careful description of the electron-phonon 
scattering processes, especially in proximity of the Fermi 
surface. This requirement translates into the necessity 
of evaluating electronic eigenvalues e k , phonon frequen- 
cies Wqi/, and electron-phonon matrix elements (7 kk v for 
a very large set of electron and phonon wavevectors in 
the Brillouin zone, of the order of tens of thousands. 

While it is practically impossible to evaluate so many 
electron-phonon matrix elements directly using standard 
density-functional calculations, it is possible to perform 
an optimal ab-initio interpolation of the matrix elements 
by exploiting localization in real space. The key idea is to 
first evaluate a small number of electron-phonon matrix 
elements in the maximally-localized Wannier representa- 
tion, and then perform a generalized Fourier interpola- 
tion into the momentum space, i.e. into the Bloch repre- 
sentation. The relation between the matrix elements in 
the Wannier representation <?rr and those in the Bloch 
representation t/ kk / is: 



(41) 



R,R' 



where N is the size of the discrete Brillouin-zone mesh, 
t/ k is a band-mixing matrix which maps electron Bloch 
bands into Wannier functions, u q is a branch-mixing ma- 
trix which maps phonon branches into individual atomic 



G 



displacements, q = k'— k, and R, R' are vectors of the di- 
rect lattice. In Eq. (|4ip the band and branch indices are 
absorbed in k, k' and in R, R'. More detailed expressions 
for implementation purposes can be found in Refs . 12711321 . 
Once obtained the matrix elements in the Wannier rep- 
resentation, the evaluation of Eq. (|4T|) for any pairs of 
initial and final electron wavevectors is inexpensive since 
it involves only very small matrix multiplications. 

The matrix elements in the Wannier representation are 
computed by hrst calculating the corresponding elements 
in the Bloch representation on a coarse Brillouin zone 
mesh using density-functional perturbation theory 4 ^ and 
then transforming into the maximally localized Wannier 
representatio n 47 ! 48 using the inverse relation of Eq. (|4T|) . 
All our density-functional and density-functional per- 
turbation theory calculations are performed using the 
Quantum ESPRESSO package^ and maximally-localized 
Wannier functions are determined using the Wannier90 
program.— The subsequent electron-phonon interpola- 
tion is performed using the EPW program,— which ex- 
tracts and processes information from both Quantum 
ESPRESSO and Wannier90. Further details on the no- 
tion of Wannier interpolation and its use in the study of 
electron-phonon interactions can be found in Ref. [28| and 
Rcfs. [27II32I respectively. 

Even when using electron-phonon Wannier interpola- 
tion the computational workload can become quite sub- 
stantial when one evaluates hundreds of thousands of ma- 
trix elements. In order to reduce the computational load 
we exploit the crystal symmetries and only evaluate the 
gap function A(k, iuJ„) and the renormalization function 
Z(k,ioj n ) in the irreducible wedge of the Brillouin zone. 
On the other hand, the sums over k' in Eqs. (f2~T|) . (|22 j) 
and Eqs. (|26I) . ([2"T)) are performed over the entire Brillouin 
zone. The meshes of wavevectors k and q= k'— k are cho- 
sen to be uniform and commensurate, in such a way that 
the grid of electron wavevectors in the final state k' maps 
into the grid of the initial wavevectors k. Since the con- 
tributions to the superconducting gap arising from elec- 
tronic states away from the Fermi energy are essentially 
negligible, the matrix elements of Eq. (j4"Tj) are evaluated 
only for electronic states such that and ek' are near the 
Fermi energy. Numerical convergence can be achieved 
typically by restricting the sums in Eqs. (j2"Tj) . (|22p and 
Eqs. (f2l) |) . (f2~T)) to an energy window around the Fermi 
level of width corresponding to 4-10 times the character- 
istic phonon frequency. 



B. Self-consistent solution of the nonlinear system 
and analytic continuation 

In order to solve iteratively the Eliashberg equations on 
the imaginary axis Eqs. (|2"TT) . (|2"2"|) . we start from an initial 
guess Ao(iui n ) for the superconducting gap. The starting 
guess A (io; n ) is chosen to be a step function vanishing 
for iu n > 2w?Jf x , Wp} 1 ax being the largest phonon energy in 
the system. The magnitude of Aq (iu! n ) is estimated from 



the BCS formula^ at zero temperature 2Ao(iw„)/T c = 
3.52, with T c given by Allen-Dynes equation^ 

Our experience shows that the convergence of the iter- 
ative self-consistent solution is significantly improved by 
using the Broydcn mixing scheme commonly employed in 
standard density-functional calculations! 52 ' 53 For the test 
cases considered in Sec. lIVI below we find that around 15- 
20 iterations are sufficient to achieve convergence when- 
ever T < 0.8T C . The number of iterations increases to 
40-60 for temperatures between 0.8-0.95T c , and may ex- 
ceed 100 for T > 0.95T C . In order to accelerate the con- 
vergence we use the gap functions calculated at a given 
temperature as seeds for the iterations at the next tem- 
perature. An alternative strategy for solving the equa- 
tions when T ~ T c would be to use the linearized form of 
the Eliashberg equations and determine the critical tem- 
perature by solving an eigenvalue problem,— 41 however 
we did not explore this possibility. 

In order to determine the superconducting gap along 
the real energy axis we consider two possibilities. The 
first possibility is to perform an approximate analytic 
continuation using Pade functions! 42 ' 43 ' 45 This procedure 
works well if the Pade functions are constructed using the 
Matsubara frequencies on the imaginary axis. The sec- 
ond possibility consists of performing the exact analytic 
continuation of the imaginary solution to the real energy 
axis as described in Sec. IIII Bl Since this latter approach 
is computationally very demanding, we speed up the con- 
vergence of the iterative analytic continuation by using 
the approximate Pade continuation as an initial guess. 



IV. APPLICATIONS 

In order to validate the computational methodology 
developed within EPW, we investigate two prototypical 
systems: the nearly-isotropic lead (Pb) superconductor, 
and the anisotropic magnesium diboride (MgB2) super- 
conductor. 



A. Computational details 

The calculations are performed within the local density 
approximation (LDA) to density-functional theor y 54 ' 55 
using Quantum ESPRESSO 42 The valence electronic wave- 
functions are expanded in planewaves basis sets with ki- 
netic energy cutoff of 80 Ry and 60 Ry for Pb and MgB 2 , 
respectively. The core-valence interaction is taken into 
account by using norm-conserving pscudopotcntialsi 56 ' 57 
For Pb we consider four valence electrons and a scalar- 
rclativistic pscudopotential. In order to facilitate the 
comparison with previous theoretical studies we use the 
LDA theoretical lattice parameters for Pb (a = 4.778 A) 
and the experimental lattice parameters for MgB2 (a = 
3.083 A and c/a = 1.142). The charge density is com- 
puted using T-centered Brillouin-zone mesh with 16 3 
and 24 3 k-points for Pb and MgB2, respectively, and a 
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FIG. 1: (Color online) Calculated isotropic Eliashberg spec- 
tral function o?F (black solid line) of Pb, and cumulative 
contribution to the electron-phonon coupling strength A (red 
dashed line). The top of the red dashed curve corresponds to 
A=1.24. 



Methfcsscl-Paxton smearing 5 ^ of 0.02 Ry. The dynami- 
cal matrices and the linear variation of the self-consistent 
potential arc calculated within density-functional pertur- 
bation theory^ 6 - on the irreducible set of a regular 8 3 (Pb) 
and 6 3 (MgB 2 ) q-point meshes. The electronic wavefunc- 
tions required for the Wannier interpolation within EPW 
are calculated on uniform and T-centered k-point meshes 
of sizes 8 3 and 6 3 for Pb and MgB 2 . 

In the case of Pb four Wannier functions are used to 
describe the electronic structure near the Fermi level. 
These states are sp 3 -like functions localized along each 
one of the Pb-Pb bonds, with a spatial spread of 2.40 A. 
In the case of MgB 2 we consider five Wannier functions 
in order to describe the band structure around the Fermi 
level. Two functions are p z -\iks states and are associated 
with the B atoms, and three functions arc cr-likc states 
localized in the middle of B-B bonds. The spatial spread 
of the MLWFs in MgB 2 are 2.02 A (p z ) and 1.16 A (a). 

In order to solve the Eliashberg equations we evalu- 
ate electron energies, phonon frequencies, and electron- 
phonon matrix elements on fine grids using the method 
of Ref. H3. The fine grids contain (40 3 ,40 3 ) k- and q- 
points for Pb (random grids), and (60 3 ,30 3 ) points for 
MgB 2 (uniform T-centered grids). Such an extremely 
fine sampling of the Brillouin zone is found to be cru- 
cial for the convergence of the superconducting energy 
gap in the fully anisotropic case. The frequency cut- 
off uj c in Eqs. ([21 ]) . ([22 ]) and Eqs. ([3"5" ]) .(|3"6 ]) is set to ten 
times the maximum phonon frequency of the system: 
uj c = lOWpJf*. The calculations are performed using 
smearing parameters in the Dirac delta functions cor- 
responding to 100 meV and 50 meV for electrons and 
phonons, respectively. 



B. Lead 

Bulk lead is the best known example of a strong- 
coupling superconductor, exhibiting a superconducting 
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FIG. 2: (Color online) Calculated energy-dependent super- 
conducting gap of Pb at T—0.3 K. The gap is obtained by 
solving the isotropic Eliashberg equations with ^ = 0.1. (a) 
Superconducting gap along the imaginary energy axis (black 
solid line), (b) Superconducting gap along the real energy 
axis. We show both the solutions obtained from the approxi- 
mate analytic continuation using Pade functions (black solid 
line and red dashed line), and the solutions obtained using 
the iterative analytic continuation (green dash-dotted line and 
blue dotted line). 



transition temperature T c = 7.2 Ki^ 9 - Although Pb is 
known to be a two-band superconductor, the supercon- 
ducting gap function is only very weakly anisotropic ^2r— 
therefore for the sake of testing our method we use the 
isotropic approximation to the Migdal-Eliashbcrg formal- 
ism described in Sec. Ill Dl 

Figure Q] shows the calculated Eliashberg spectral func- 
tion a 2 F(ui) and the corresponding electron-phonon cou- 
pling parameter A. We find an overall good agreement 
with experimental results^, although wc observe a small 
(~0.5 meV) but non-negligible blueshift of the two peaks 
in the Eliashberg function. This blueshift is well un- 
derstood now and arises from the overestimation of the 
phonon frequencies in the absence of spin-orbit coupling 
in our calculation ! 64 ' 65 Our calculated electron-phonon 
coupling A=1.24 lies in between the values reported in 
previous theoretical studies ; 39 ' 62 ' 63 although it is some- 
what smaller than the value 1.55 obtained from tunneling 
measurements^ 9 - owing to the neglect of spin-orbit cou- 
pling. 

Figure [2] shows the solutions of the isotropic Eliashberg 
equations Eqs. (|55 ]) .([5rJ ]l for /i* = 0.1 and T=0.3 K. Along 
the imaginary axis the superconducting gap function is 
purely real and displays a frequency dependence similar 
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FIG. 3: (Color online) (a) Calculated quasiparticle density 
of states of Pb at T—0.3 K (black solid lines). The supercon- 
ducting gap is obtained from Fig. [2] (b) Same quantity as in 
(a), magnified in order to show the structure which is used in 
tunneling experiments for extracting the Eliashberg spectral 
function. 



to standard plasmon-pole models [Fig.[2ja)]. The contin- 
uation of the calculated superconducting gap to the real 
energy axis is shown in Fig.[3Jb). We see that the approx- 
imate analytic continuation using Pade functions and the 
exact iterative analytic continuation yield very similar 
results. As expected the approximate Pade continuation 
misses some of the fine features which are instead ob- 
served in the exact iterative continuation. Our calculated 
superconducting gap is in very good agreement with solu- 
tions of the Eliashberg equations obtained directly on the 
real energy axis<£2 In Fig. [3Jb) we see that a two-peak 
structure emerges both in the real part and the imagi- 
nary part of A (a;). These two peaks occur on the scale 
of the phonon energies and correlate with those observed 
in the Eliashberg spectral function of Fig. [T}£ A detailed 
analysis shows that the peaks in the real part of the gap 
function are blueshifted by approximately Ao = A(cj = 0) 
with respect to the corresponding peaks in a 2 F(u). 

Figure 03a) shows the normalized quasiparticle density 
of unoccupied states obtained from Eq. (|40"]) using the 
gap function of Fig. [2] As expected the strong van Hove 
singularity marks the leading edge Ao of the supercon- 
ducting gap. The fine structure of the density of states 
around the van Hove singularity [Fig. \5jlb)} is the direct 
signature of the electron-phonon physics and is precisely 
the basis for direct measurements of the Eliashberg func- 
tion using tunneling spectroscopy. 

FigurcUshows the superconducting gap function at the 




FIG. 4: (Color online) (a) Calculated superconducting gap 
of Pb at the Fermi level as a function of temperature (disks). 
The Coulomb parameter is set to ii* = 0.1. (b) Calculated 
superconducting gap of Pb at the Fermi level for T=0 K as a 
function of the Coulomb parameter fj,* (disks), (c) Calculated 
critical temperature as a function of the Coulomb parameter 
fil (disks). In all panels the solid lines are guides to the eye. 



Fermi level as a function of temperature, calculated for 
fi* =0.1. The leading edge of the gap at T=0 K is found 
to be Ao=1.24 meV, in good agreement with tunneling 
measurements yielding 1.33 meV.— The superconduct- 
ing T c is identified as the temperature at which the gap 
vanishes. From Fig.[4|a) we find T c =6.8 K, in very good 
agreement with previous theoretical studies,—"^ and only 
slightly lower than the experimental datum of 7.2 K. For 
completeness in Fig.@|b) and (c) we also explore the sen- 
sitivity of the calculated gap and critical temperature to 
the choice of the Coulomb parameter fx* . As expected, a 
reduction of the effective Coulomb interaction results in 
an increase of both Aq and T c . 



C. Magnesium diboride 

Within the class of phonon-mediated superconductors 
MgB2 holds the record of the highest critical tempera- 
ture, with T c =39 Kj££ After a decade of intense exper- 
imental and theoretical investigations since its discov- 
ery, it is now understood that MgB2 is an anisotropic 
two-gap electron-phonon superconductori 15 ' 17 i 67 ~ — The 
anisotropy of the superconducting gap is a consequence 
of the multi-sheet Fermi surface of MgB 2 , consisting of 
two hole-like coaxial cylinders arising from the a bonding 
bands, and two hole-like tubular networks arising from 
the 7T bondin g a nd antibonding bands (see for example 
Fig. 3 of Ref.ES). 

Figure 03a) shows our calculated isotropic Eliash- 
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FIG. 5: (Color online) (a) Calculated isotropic Eliashberg 
spectral function o?F of MgB2 (black solid line) and cumu- 
lative contribution to the electron-phonon coupling strength 
A (red dashed line), (b) Distribution of the electron-phonon 
coupling strenght Ak for MgB2 (black solid line). 



berg spectral function a 2 F(ui) and electron-phonon cou- 
pling strength A for MgB2. a 2 F(ui) displays a large 
dominant peak centered around 60 meV and a sec- 
ond weaker peak centered around 86 meV. The corre- 
sponding isotropic electron-phonon coupling strength is 
A=0.748. In order to quantify the different contributions 
to the coupling strength associated with the a sheets and 
from the it sheets of the Fermi surface we evaluate a 
band- and wavevector- dependent electron-phonon cou- 
pling strength defined by Ak = J2k< Wk'A(k,k',n = 0). 
Figure [SJb) shows that the calculated Ak cluster into two 
separate ranges. The lower range Ak=0.35-0.50 corre- 
sponds to the coupling of the n Fermi surface sheets, and 
the higher range 0.95-1.35 corresponds to the coupling of 
the a sheets. The wider range of Ak in the a sheets re- 
flects a more pronounced anisotropy with respect to the ir 
sheets. The structure of a 2 F(uj) and the calculated cou- 
pling strength are in good agreement with previous cal- 
culations J^ii2iZ2r 1 22 j n p ar ticular our results are in very 
good agreement with those reported in Ref. [7l| where a 
related interpolation scheme was used. For completeness 
we show in Table U the sensitivity of the calculated aver- 
age coupling strength A on the underlying Brillouin-zone 
grids, and we compare with previous first-principles cal- 
culations. 

Figure[6ja) shows the anisotropic superconducting gap 
function A(k, uj) of MgB2 at T=10 K calculated along 
the imaginary axis using the anisotropic Eliashberg equa- 
tions Eqs. 421]), and ^*=0.16. Figure [IJb) shows the 



real part of the superconducting gap function along the 
real energy axis, as obtained from the imaginary axis 
solutions of Fig EJa) via the approximate analytic con- 
tinuation using Padc functions. The two-gap nature of 
MgB2 emerges in a completely natural way from our im- 
plementation. Indeed for each energy two distinct sets 
of superconducting gaps can be identified and associated 
with the a and the it sheets of the Fermi surface. The two 
gaps are both anisotropic, and the corresponding Fcrmi- 
surface averages are A T = 1.8 meV and A a = 8.5 meV, 
respectively. For comparison the experimental values for 
the gaps lie in the range 2.3-2.8 meV for 7t-band, and 7.0- 
7.1 meV for a-bandi^r— As in the case of Pb, the struc- 
ture that can be observed both in the a and ir supercon- 
ducting gaps reflect the peaks occuring in the Eliashberg 
spectral function of Fig. [5] 

Figure [TJa) shows the calculated leading edges Ak of 
the superconducting gaps as a function of temperature. 
Both the 7T and a gaps vanish at the critical tempera- 
ture T c =50 K. The corresponding quasiparticlc density 
of states, presented in Fig. [TJb), clearly shows the two- 
gaps structure of MgB2, in agreement with experiment ^ 

Our calculated critical temperature is larger than the 
experimentally measured T c of 39 however it is in 
very good agreement with previous first-principles cal- 
culations based on the ME formalism-^ or the SCDFT 
formalism^ At this time it is still unclear whether the 
overestimation of the experimental critical temperature 
is due to possible anharmonic effects^ or to the use 
of an isotropic Coulomb parameter^ In fact, using the 
anisotropic Eliashberg formalism and a Coulomb param- 
eter /x* = 0.12, the authors of Ref. [I?] find that the cal- 
culated T c decreases from 55 K to 39 K if phonon anhar- 



Reference 


k-mesh 


q-mesh 


A 


Bohnen et al. (Ref. 70) 


36 3 


6 a 


0.73 


Choi et al. (Ref. .17) 


12xl8 2 


12 x 18 2 


0.73 


Floris et al. (Ref. _15) 


24 3 


20 3 


0.71 


Eiguren et al. (Ref. 71) 


40 3 


40 3 


0.776 


Calandra et al. (Ref. 72) 


80 3 


20 3 


0.741 


This work 


40 3 


20 3 (40 3 ) 


0.735 




80 3 


20 3 (40 3 ) 


0.739 




60 3 


30 3 (60 3 ) 


0.748 




30 3 


30 3 


0.782 




50 3 


50 3 


0.744 




64000 


8000 


0.757 




216000 


27000 


0.726 



TABLE I: Electron-phonon coupling strength A of MgB2 cal- 
culated using various meshes of k- and q-points in the Bril- 
louin zone. The numbers in the brackets correspond to a sec- 
ond choice of q-mesh while keeping the k-mesh unchanged. 
The two bottom rows correspond to uniform random distri- 
butions of k- and q-points. 
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FIG. 6: (Color online) Calculated energy-dependent super- 
conducting gap of MgB2 at T=10 K. The gap is obtained 
by solving the fully anisotropic Eliashberg equations with 
He = 0.16. (a) Superconducting gap along the imaginary 
energy axis (black dots), (b) Superconducting gap along the 
real energy axis, obtained from the approximate analytic con- 
tinuation using Pade functions (black dots) . Only a represen- 
tative sample of 10 data points out of entire set of calculated 
gaps (10 7 points) is shown for clarity. 



monicity is taken into account. On the other hand, using 
the SCDFT formalism the authors of Ref. [TH calculate a 
critical temperature T c =22 K when using the complete 
wavevector-dependent superconducting gap. When em- 
ploying a band-averaged superconducting gap and var- 
ious levels of approximations for the Coulomb interac- 
tion, the same authors find critical temperatures in the 
range 30-50 K.— A similar sensitivity of the calculated 
T c to fine details of the calculations are reported in other 
studies based on a two-bands approximation to the ME 
formalism.i^ The origin of the discrepancy between 
first-principles calculations of the critical temperature of 
MgB2 and experiment clearly deserves further investi- 
gation, however this is beyond the scope of the present 
manuscript. 

The use of electron-phonon Wannicr interpolation al- 
lows us to investigate the sensitivity of the superconduct- 
ing gap to the electron and phonon meshes used for the 
calculations. Figure [5] shows the energy distribution of 
the cr-gap at T=30 K for eight sets of electron and phonon 
meshes. If we compare the gaps shown in Fig. [3] and the 
average coupling strengths reported in TableUwe see that 
obtaining converged results for Ak is considerably more 
challenging than for A. For example, when using the same 
k-points mesh (80 3 ) and different q-points meshes (20 3 



FIG. 7: (Color online) (a) Calculated anisotropic supercon- 
ducting gaps of MgB2 on the Fermi surface as a function of 
temperature. The Coulomb potential is set to [i* — 0.16. (b) 
Corresponding quasiparticle density of states for a few repre- 
sentative temperatures (15 K black solid line, 30 K red dashed 
line, 45 K blue dash-doted line). 



and 40 3 ), the spread of the superconducting gap distri- 
bution changes from ~2.5 meV to ~1.5 meV, while the 
average coupling strength A is essentially unaffected. The 
same observation applies when we compare results for the 
same q-points mesh (40 3 ) but different k-points meshes 
(40 3 and 80 3 ). These differences highlight the difficulty of 
describing anisotropic quantities, and point out the need 
of having a very dense sampling of the Brillouin zone not 
only for the electrons but also for the phonons. For this 
reason the combination of the Eliashberg formalism with 
electron-phonon Wannicr interpolation demonstrated in 
this work provides an ideal computational tool for inves- 
tigating anisotropic superconductors. 



V. CONCLUSIONS 

In summary we developed a computational method 
which combines the anisotropic Migdal-Eliashbcrg for- 
malism with electron-phonon interpolation based on 
maximally-localized Wannicr functions (EPW). Our new 
method allows us to calculate the momentum- and band- 
resolved superconducting gap both effectively and accu- 
rately using a very fine description of electron-phonon 
scattering processes on the Fermi surface. In order to 
demonstrate our methodology we reported a compre- 
hensive set of tests on two representative superconduc- 
tors, namely Pb and MgB2, and validated our approach 
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FIG. 8: (Color online) Calculated energy distribution of the 
superconducting gaps on the a sheets of the Fermi surface of 
MgB2 at T=30 K (black lines). The gaps are obatined for 
various k- and q-points meshes, as indicated by the labels 
on the horizontal axis. The Coulomb parameter is set to 
ill = 0.16. 

against previous first-principles calculations as well as ex- 
periment. We discussed the performance of two analytic 
continuation methods for obtaining the superconducting 
gap on the real energy axis, and we investigated the sen- 
sitivity of the calculated gaps on the underlying choice 
of the Brillouin-zone grids for electrons and phonons. 

In order to set a road map for first-principles studies of 
superconductors we now discuss the key approximations 
involved in the present approach and suggest possible 
avenues for future developments. The main approxima- 
tions leading to Eqs. (|21[) . ([22]) arc as follows: (i) vertex 
corrections in the diagrammatic expansion of the electron 
self-energy are neglected following Migdal's theorem, (ii) 
the self-energy is assumed to be diagonal in band indices, 
(iii) the Eliashberg equations are restricted to bands near 
the Fermi level, (iv) the density of states near the Fermi 
level is assumed to be constant, and (v) the Coulomb 
interaction is described by an empirical parameter i 2 ' 35 i 36 

Regarding approximation (i), it is well known that 
Migdal's theorem can break down for large electron- 
phonon coupling or in the presence of Fermi-surface nest- 
mg££5- In this area the availability of electron-phonon 
matrix elements at a very small computational cost, as 
provided by our method, could be used as a starting 
point to explore the effects of vertex corrections beyond 
Migdal's approximation. For example the evaluation of 
the first non-crossing diagram should not constitute a 
major challenge, at least in the normal state. This may 
help assessing the numerical error introduced by Migdal's 
approximation. 

Going beyond approximation (ii) is computationally 



challenging since the Green's function and the self-energy 
should be treated as matrices in the band indices i 76 i 77 
This step would increase substantially the complexity of 
the formalism since the inversion of the Dyson's equa- 
tion would require matrix operations. Given the small 
energies associated with the superconducting gap it is 
reasonable to expect that off-diagonal terms should not 
play an important role in simple cases. However in the 
presence of degeneracies, and in particular in Jahn- Teller 
systems, the correct description of these terms may prove 
critical^ 

It should be possible, at least in principle, to remove 
approximation (iii) by including bands away from the 
Fermi level in the calculations. Here the difficulty is 
purely on the computational side, and for simple systems 
this should be doable. However relaxing this constraint 
would also introduce one additional equation [Eq. (|16|) ] 
in order to calculate the correct quasiparticle shifts and 
impose the conservation of charge in the system. 

The assumption (iv) of a constant density of states 
near the Fermi level obviously breaks down for materi- 
als exhibiting structure in the density of states on the 
scale of the phonon energy. In theses cases it has been 
shown that the energy dependence of density of states 
can be retained within the isotropic approximation to 
the Eliashberg equationsi 2 ' 35 i 79 We are planning to in- 
clude this possibility in our methodology in the future. 

Lastly, approximation (v) means that in the present 
implementation the Coulomb repulsion between electrons 
remains described at the empirical level as an adjustable 
parameter. In order to remove this limitation our first 
step will be to evaluate the Coulomb parameter using the 
dielectric matrix in the random-phase approximation^ 
In the longer term it would be desirable to introduce in- 
side Eq. ([212]) matrix elements of the screened Coulomb in- 
teraction calculated using the Sternheimer-GW method 
of Ref. iH. 

We hope that the method reported here will prove use- 
ful to the superconductivity community as a robust and 
rigorous procedure for shedding light on existing super- 
conductors and possibly for predicting new superconduc- 
tors yet to be discovered. 
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